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Abstract. We report on numerical simulations of the two-dimensional Blume-Capel ferromagnet embedded 
in the triangular lattice. The model is studied in both its first- and second-order phase transition regime for 
several values of the crystal field via a sophisticated two-stage numerical strategy using the Wang-Landau 
algorithm. Using classical finite-size scaling techniques we estimate with high accuracy phase-transition 
temperatures, thermal, and magnetic critical exponents and we give an approximation of the phase diagram 
of the model. 

PACS. PACS. 05.50+q Lattice theory and statistics (Ising, Potts, etc.) - 64.60.Fr Equilibrium properties 
near critical points, critical exponents - 75.10.Hk Classical spin models 

1 Introduction malization groups [6], high- and low-temperature series 

calculations [7j , a phenomenological finite-size scaling anal- 
The Blume-Capel model consisting of a spin-one Ising ysis using a strip geometry [819] , and, finally, the most 
Hamiltonian with a single-ion uniaxial crystal field anisotropyrffjajit numerica l approaches via the Wang-Landau algo- 
is one of the most studied models in the communities rithm [10111112] . The phase diagram of the model consists 

of a segment of continuous Ising-like transitions at high 
temperatures and low values of the crystal field which 
ends at a tricritical point, where it is joined with a sec- 
ond segment of first-order transitions between (A t ,T t ) and 
(A = z J/2, T = 0), where the subscript t refers to the tri- 
critical point and z denotes the coordination number of 
the considered lattice. 

In the present paper we are interested in the critical 
properties of the Blume-Capel model embedded in the 
triangular lattice {z — 6). For this particular case, Ma- 
han and Girvin [T3] were the first to apply position-space 
renormalization-group methods to solve the model. These 
authors estimated the critical frontier of the model with 
the location of the tricritical point (A t ,T t ) at (2.686, 1.493). 
Many years later, Du et al. 14 performed a further sophis- 
ticated analytical calculation using an expanded Bethe- 
Peierls approximation to find respectively the estimates 
(A t ,T t ) = (2.818,1.598). Here, we employ MC simula- 
tions to investigate several aspects of the phase diagram 
of the model and refine mean-field-type estimates, which, 
as is well known, suffer from large errors. Q Our exten- 
sive simulations follow a sophisticated numerical scheme, 
as described in the following Section, that enabled us to 
perform a detailed finite-size scaling analysis of the critical 
properties of the model in both its first- and second-order 
phase transition regimes. These results, together with an 



of Statistical Mechanics and Condensed Matter Physics. 
This is not only because of the relative simplicity with 
which approximate calculations for this model can be car- 
ried out and tested, as well as the fundamental theoreti- 
cal interest arising from the richness of its phase diagram, 
but also because versions and extensions of the model can 
be applied for the description of many different physical 
structures, some of them being multi-component fluids, 
ternary alloys, and 3 He - 4 He mixtures [3]. Noteworthy, 
latest applications of the Blume-Capel model include an- 
alyzes of ferrimagnets, as recently discussed in a thorough 
contribution by Selke and Oitmaa [3]. 

The Blume-Capel model is described by the Hamilto- 
nian 

H = -J ]T s iSj + Aj2sl (1) 

<ij> i 



where the spin variables Si take on the values —1,0, or 
+1, < ij > indicates summation over all nearest-neighbor 
pairs of sites, and J > is the ferromagnetic exchange in- 
teraction (here we set J = 1 and ks — 1 to fix the temper- 
ature scale). The parameter A is known as the crystal- field 
coupling that controls the density of vacancies (si — 0). 
For A — > — oo vacancies are suppressed and the model 
maps onto the Ising model. 

As it is well known, the model of equation ([1} has 
been analyzed, besides the original mean- field theory [112] . 
by a variety of approximations and numerical approaches, 1 Note here that the estimation of Du et al. for the corre- 
mostly on the square lattice. These include the real space sponding tricritical crystal field value of the square lattice is 
renormalization group, Monte Carlo (MC) simulations and A t = 1.906, whereas the most accurate MC estimation in the 
MC renormalization-group calculations [5], e-expansion renorliterature is A t = 1.96581... by Qian et al. [15] . 
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approximation of the phase diagram of the model and the 
estimation of the tricritical value of A t via a new scheme 
are presented in Section [3] This contribution is ended in 
Section 0] where a brief summary of our conclusions is 
given. 

2 Outline of the Numerical Approach 

In the last few years we have used an entropic sampling 
implementation of the Wang-Landau algorithm [TB] to 
study some simple [17], but also some more complex sys- 
tems [18119120] . One basic ingredient of this implementa- 
tion is a suitable restriction of the energy subspace for 
the implementation of the Wang-Landau algorithm. This 
was originally termed as the critical minimum energy sub- 
space restriction |17] and it can be carried out in many 
alternative ways, the simplest being that of observing the 
finite-size behavior of the tails of the energy probability 
density function of the system [IT] . 

Complications that may arise in complex systems, i.e. 
random systems or systems showing a first-order phase 
transition, can be easily accounted for by various simple 
modifications that take into account possible oscillations 
in the energy probability density function and expected 
sample-to-sample fluctuations of individual realizations. 
In our recent papers [12119120) . we have presented de- 
tails of various sophisticated routes for the identification of 
the appropriate energy subspace (E\,E2) for the entropic 
sampling of each realization. In estimating the appropri- 
ate subspace from a chosen pseudocritical temperature 
one should be careful to account for the shift behavior 
of other important pseudocritical temperatures and ex- 
tend the subspace appropriately from both low- and high- 
energy sides in order to achieve an accurate estimation of 
all finite-size anomalies. Of course, taking the union of the 
corresponding subspaces, insures accuracy for the temper- 
ature region of all studied pseudocritical temperatures. 

The up to date version of our implementation uses a 
combination of several stages of the Wang-Landau pro- 
cess. First, we carry out a starting (or preliminary) multi- 
range (multi-R) stage, in a very wide energy subspace. 
This preliminary stage is performed up to a certain level 
of the Wang-Landau random walk. The Wang-Landau re- 
finement is G(E) — > f-G(E), where G{E) is the density of 
states (DOS) and we follow the usual modification factor 
adjustment fj+i = \J~fj and /i = e. The preliminary stage 
may consist of the levels : j = = 18 and to im- 

prove accuracy the process may be repeated several times. 
However, in repeating the preliminary process and in or- 
der to be efficient, we use only the levels j = 13, . . . , 18 
after the first attempt, using as starting DOS the one ob- 
tained in the first random walk at the level j — 12. From 
our experience, this practice is almost equivalent to sim- 
ulating the same number of independent Wang-Landau 
random walks. Also in our recent studies we have found 
out that is much more efficient and accurate to loosen 
up the originally applied very strict flatness criteria [17] . 
Thus, a variable flatness process starting at the first lev- 
els with a very loose flatness criteria and assuming at the 



level j = 18 the original strict flatness criteria is nowa- 
days used. After the above described preliminary multi-R 
stage, in the wide energy subspace, one can proceed in 
a safe identification of the appropriate energy subspace 
using one or more alternatives outlined in reference [17] . 

The process continues in two further stages (two-stage 
process), using now mainly high iteration levels, where 
the modification factor is very close to unity and there is 
not any significant violation of the detailed balance con- 
dition during the Wang-Landau process. These two stages 
are suitable for the accumulation of histogram data (for 
instance energy- magnetization histograms), which can be 
used for an accurate entropic calculation of non-thermal 
thermodynamic parameters, such as the order parameter 
and its susceptibility [T7]. In the first (high-level) stage, we 
follow again a repeated several times (typically ~ 5 — 10) 
multi-R Wang-Landau approach, carried out now only in 
the restricted energy subspace. The Wang-Landau levels 
may be now chosen as j = 18, 19, 20 and as an appropriate 
starting DOS for the corresponding starting level the aver- 
age DOS of the preliminary stage at the starting level may 
be used. Finally, the second (high-level) stage is applied in 
the refinement Wang-Landau levels j = ji, ■ ■ ■ ,ji + 3 (typ- 
ically ji = 21), where we usually test both an one-range 
(one-R) or a multi-R approach with large energy inter- 
vals. In the case of the one-R approach we have found 
very convenient and in most cases more accurate to fol- 
low the Belardinelli and Pereyra [21] adjustment of the 
Wang-Landau modification factor according to the rule 
In / ~ i -1 . Finally, it should be also noted that by apply- 
ing in our scheme a separate accumulation of histogram 
data in the starting multi-R stage (in the wide energy sub- 
space) offers the opportunity to inspect the behavior of all 
basic thermodynamic functions in an also wide tempera- 
ture range and not only in the neighborhood of the finite- 
size anomalies. The approximation outside the dominant 
energy subspace is not of the same accuracy with that 
of the restricted dominant energy subspace but is good 
enough for the observation of the general behavior and 
provides also a route of inspecting the degree of approxi- 
mation. 



In the present work, the above described numerical 
approach was used to estimate the properties of the tri- 
angular Blume-Capel model for lattice sizes in the range 
L E {20 — 100} for all values of the crystal field con- 
sidered. For each pair {L, A), 50 independent runs were 
performed. We close this outline of our numerical scheme 
with some comments concerning statistical errors. Even 
for the larger lattice size studied here (L = 100), and de- 
pending on the thermodynamic parameter, the statistical 
errors of the Wang-Landau method were found to be of 
reasonable magnitude and in some cases to be of the order 
of the symbol sizes, or even smaller. Thus, the error bars 
shown in all our figures in the following Section have been 
estimated as standard deviations from the ensemble of the 
50 independent runs. 
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Fig. 1. Energy probability density functions for a lattice size 
L — 60 and several values of the crystal field, as indicated. 
For the cases of a double-peaked structure we have used the 
temperature where the two peaks are of equal height, whereas 
in the other single-peaked cases, we plot the densities at the 
corresponding pseudocritical temperatures of the specific heat. 

3 Results and Analysis 

We present in this Section our numerical results and analy- 
sis for the triangular Blume-Capel model. For convenience 
we separate our discussion in three parts. The first (sub- 
section 13. 1[) refers to the second-order phase transition 
regime of the model, the following (subsection I3.2j) to the 
first-order phase transition regime, and finally in the last 
subsection l3.3l we give our approximation on the phase di- 
agram of the model together with a novel estimation ap- 
proach of the tricritical crystal field value A t . Preliminary 
runs, as shown in figures [1] and [51 indicate that the loca- 
tion of the tricritical point of the model is in the regime 
A > 2.915. More specifically, in figurc[T]we plot the energy 
probability density function P(e), where e = H/L d and 
d the lattice dimensionality (d = 2 in the present study), 
for a lattice size L = 60 for several values of the crystal 
field in the regime A = 2.9 — 2.95. The double-peaked 
structure in the energy probability density function, sig- 
naling a first-order phase transition |22l23j , appears after 
the value A = 2.915. Respectively, figure [5] illustrates the 
corresponding specific heats C (main panel) and magnetic 
susceptibilities \ (inset) as a function of the temperature 
for the same values of the crystal field of figure [TJ Again 
the sharp peak in both quantities, characteristic of a first- 
order phase transition, is observed in the same Z\-regime. 

Thus, using the above information, we have chosen to 
simulate the values A e {0,0.5,1,1.5,2,2.5,2.75} in the 
second-order phase transition regime of the model and the 
value A = 2.95 in the first-order phase transition regime, 
respectively. 



3.1 Second-order phase transition regime 

The triangular Blume-Capel model at the crystal field 
value A < A t , undergoes a second-order transition be- 



Fig. 2. Corresponding specific heats (main panel) and mag- 
netic susceptibilities (inset) as a function of temperature. This 
figure is complementary to figure [T] 



tween the ferromagnetic and paramagnetic phases, ex- 
pected to be in the universality class of the simple d = 2 
Ising model. In the following, we present the finite-size 
scaling analysis of our numerical data for this case, to 
verify this expectation. In figure[3]we give an example, for 
the case A — 1, of the shift behavior of the pseudocritical 
temperatures corresponding to the peaks of the following 
six quantities: specific heat C, magnetic susceptibility x, 
derivative of the absolute order parameter with respect to 
the inverse temperature K = 1/T [21] 



d(\M\) 
dK 



= (\M\H) (\M\)(H), 



(2) 



and logarithmic derivatives of the first-, second-, and fourth- 
order powers of the order parameter with respect to the 
inverse temperature K = 1/T [2~4] 



d\n(M n ) (M n H) 



(H). 



dK (M n ) 
Fitting our data to the expected power-law behavior 



(3) 



(4) 



we find the critical temperature to be T c = 2.395(2) and 
the correlation length exponent v — 1.007(19), in agree- 
ment with the value of v = 1 of the simple Ising model. 
The analysis presented in figure [3] has been performed for 
all the values of the crystal field considered in the present 
paper to estimate transition temperatures, but are not 
shown here for brevity. These values will be used in the 
sequel in the construction of the phase diagram of the 
model. 

Figure @] gives an alternative estimation of the correla- 
tion length exponent via the finite-size scaling behavior of 
the logarithmic derivatives of the order parameter, whose 
maxima are expected to scale (in a second-order phase 
transition) as ~ L}l v with the system size [24]. We chose 
to show in figure @] the case n = 2 of equation ^ . A si- 
multaneous fitting for all the values of the crystal field 
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Fig. 3. Simultaneous fitting of the form ((4| of the six pseud- 
ocritical temperatures defined in the text for the case A = 1. 



Fig. 4. Simultaneous fitting of the maxima of the second-order 
logarithmic derivative of the order parameter for A £ [0 — 2.75]. 



considered gives an excellent estimate of 1.008(11) for the 
critical exponent v of the correlation length. 

We proceed with figures [5] and [5] to estimate the mag- 
netic exponent ratios of the triangular Blume-Capel model. 
In particular in figure [5] we present the order-parameter 
data at the estimated critical temperatures as a function 
of the lattice size. The lines show a simultaneous fitting for 
the whole spectrum of Z\-values of the form M c ~ L~^l v 
that gives an estimate 0.126(2) for the critical exponent 
ratio /3/v, in very good agreement with expected value 
0.125 of the pure Ising case. Respectively, in figure [6] we 
plot the maxima of the magnetic susceptibility, again as a 
function of the lattice size. As in figure the lines show a 
simultaneous power-law fitting of the form \* ~ L 1 ^ that 
gives an estimate 1.756(9) for the magnetic exponent ra- 
tio 7/V, also in very good agreement with expected value 
1.75 of the pure Ising case. 

Closing this subsection, we deal with the most intrigu- 
ing thermodynamic quantity in the study of spin models: 
the specific heat. The specific heat is an extremely sensi- 
tive quantity and it is well known that it is, at least in some 
cases, a very hard task to identify with good accuracy its 
scaling behavior. However, in the present work, our care- 
ful numerical implementation of the Wang-Landau scheme 
and the repeated sampling gave very well numerical data, 
as can be shown in figure [JJ where we present the finite- 
size scaling behavior of the specific heat maxima. The solid 
lines show a excellent simultaneous logarithmic fitting of 
the form 



C* = d+C 2 lnL. 



(5) 



Note here that, the expected logarithmic divergence of the 
specific heat is very well obtained even from the smaller 
lattice sizes shown. 

Summarizing, the results and analysis presented in this 
subsection for the triangular Blume-Capel model at its 
second-order phase transition regime are in full agreement 
with universality arguments that place the Blume-Capel 
model for A < A t in the Ising universality class. 
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Fig. 5. Simultaneous fitting of the critical order-parameter 
data for A € [0 - 2.75]. 
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y/v = 1.756(9) 




Fig. 6. Simultaneous fitting of the maxima of the magnetic 
susceptibility for A £ [0-2.75]. 
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Fig. 7. Finite-size scaling of the specific heat maxima. The Fig. 8. The double-peaked structure of the energy probability 



lines show a simultaneous fitting of the form ([5]) for A € [0 
2.75]. Note the logarithmic scale in the horizontal axis. 



density function at T — Th for A = 2.95 and L = 60. 



3.2 First-order phase transition regime 

We now consider the case A = 2.95, for which the model 
undergoes a first-order transition between the ferromag- 
netic and paramagnetic phases. Our first attempt to elu- 
cidate the first-order transition features of the present 
model will closely follow previous analogous studies car- 
ried out on the q > 5 Potts model [25126127128] and also 
our studies of the corresponding square lattice Blume- 
Capel model [T2| and the triangular Ising model with 
nearest- and next-nearest-neighbor antiferromagnetic in- 
teractions As it is well known from the existing the- 
ories of first-order transitions, all finite-size contributions 
enter in the scaling equations in powers of the system size 
L d [3D] . This holds for the general shift behavior and also 
for the finite-size scaling behavior of the peaks of various 
energy cumulants and of the magnetic susceptibility. It 
is also well known, as mentioned above, that the double- 
peaked structure of the energy probability density func- 
tion P{e) is signaling the emergence of the expected two 
delta-peak behavior in the thermodynamic limit, for a gen- 
uine first-order phase transition |22l23j , and with increas- 
ing lattice size the barrier between the two peaks should 
steadily increase. According to the arguments of Lee and 
Kosterlitz [26] the quantity AF{L)/L d ~ l = [k B T\n (P mi 
where P max and P m in are the maximum and minimum 
energy probability density function values at the temper- 
ature T/- t where the two peaks are of equal height, should 
tend to a non-zero value. 

Figure [5] shows the pronounced double-peaked struc- 
ture of the energy probability density function of the model 
at the temperature T — Th where the two peaks are of 
equal height for L = 60. From the double-peaked energy 
probability density function one can estimate the surface 
tension S(L) — AF(L)/L and the latent heat Ae(L), 
whose values remain finite for a genuine first-order tran- 
sition. Figure [9] shows the limiting behavior of these two 
quantities and verifies the persistence of the first-order 
character of the transition at A = 2.95. The limiting val- 
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Fig. 9. Limiting behavior for A — 2.95 of the surface tension 
S(L), defined in the text, and the latent heat Ae(L), shown 
in figure [8] The solid lines are linear fittings to L — > oo. 



ues of S(L) and Ae(L) are given in the graph by extrap- 
olating at the larger lattice sizes studied. 

Figures ITUl and [TT1 illustrate that the traditionally used 
divergences in finite-size scaling of the specific heat C 
and susceptibility x follow very well a power-law behav- 
^/-ftjie form ~ L d , as expected for first-order tran- 
sitions [22 23 . Furthermore, figure IT21 demonstrates that 
the divergences corresponding to the first-, second-, and 
fourth-order logarithmic derivatives of the order parame- 
ter defined in equation Q follow also very well the same 
L d behavior, as expected. 

In this subsection we have presented a reliable analy- 
sis of the first-order transition features of the triangular 
Blume-Capel model at the value A = 2.95 of the crystal 
field and we have verified all the theoretical expectations 
of finite-size scaling in first-order phase transitions. Ad- 
ditionally, we have estimated important features of this 
transition, i.e. the surface tension and latent heat of the 
transition, using the well-established method of Lee and 
Kosterlitz [55]. 
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Fig. 10. Finite-size scaling behavior of the specific heat peaks 
at A = 2.95 giving a clear L d divergence. 
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Fig. 11. Finite-size scaling behavior of the susceptibility max- 
ima at A = 2.95. A clear L d divergence is observed as in the 
case of the specific heat. 
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Fig. 12. Simultaneous fitting of the maxima of the three loga- 
rithmic derivatives of the order parameter ((3j) for n = 1, 2, and 
4, at A = 2.95. As in figures [TOl and ITTI a clear L d divergence 
is observed. 




A= 2.925(8) 

0.0 I 1 1 ' 1 ' 1 ' 1 1 1 ■— Li 1 

0.0 0.5 1.0 1.5 2.0 2.5 3.0 3.5 
A 



Fig. 13. An approximate estimation of the tricritical value of 
the crystal field by fitting the decreasing logarithmic ampli- 
tudes Ca of the model [see equation (0] at suitable values of 
the crystal field to a power law, as shown in the figure. 

3.3 Phase diagram and tricriticality 

We close this Section with our results on the phase dia- 
gram of the model and the related issue of tricriticality. 
From figure [7J of subsection 13.11 one can observe the ex- 
pected Ising logarithmic divergence of the specific-heat 
maxima. Avoiding the value A = 0, which suffers from 
(small) crossover effects, we attempted to estimate the tri- 
critical value of the crystal field by fitting the decreasing 
logarithmic amplitudes C2 estimated in the simultaneous 
fitting of figure [7] to a suitable power law, as shown in the 
figure [T3l This may look like a questionable idea, since 
the behavior of specific heat data is the Achilles' heel of 
finite-size scaling analysis. Yet, figure H3] shows that be- 
sides the rather large errors in the logarithmic amplitudes 
C2 , one may approximately estimate the tricritical crystal- 
field value to be A t — 2.925(8), as shown in the panel of 
the figure, in excellent agreement with the previous anal- 
ysis of the energy probability density functions and the 
related thermodynamic quantities of figures [T] and [H We 
should note here that, such an idea was first suggested 
and performed in the recent papers of reference |12) for 
the cases of the pure and random-bond square Blume- 
Capel model, where it has also produced quite acceptable 
results. Furthermore, we aim here only to a qualitative 
prediction of the tricritical crystal field value. A most ac- 
curate prediction of these first-order transition features 
of the phase diagram of the model should follow different 
routes, well-established in the literature, which are, in any 
case, well beyond the scope of the present paper. 

Using now our estimates for the transition tempera- 
tures obtained throughout this paper for all the values 
of the crystal field considered we attempt to construct in 
figure [U an approximation of the phase diagram of the 
model. We have fitted the points shown in figure [F4"l using 
the following power-law ansatz 

T~(A -Af, (6) 
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Fig. 14. Approximation of the phase diagram of the triangular 
Blume-Capel model. The dotted and solid curves correspond 
respectively to first- and second-order phase transition lines 
that separate the ordered and paramagnetic phases. The black 
rhombus points the approximate location of the tricritical point 
(TP). 

where Aq in the above equation © denotes the crossing 
of the boundary to the horizontal axis (T = 0) and should 
take the value 3(= zJ/2) for the present case of the trian- 
gular lattice [z — 6). The solid and dotted curves that cor- 
respond to second- and first-order phase transitions lines 
that separate the ordered from the paramagnetic phase 
have been obtained from a fitting of the form (J6j) that 
gives the value 3.01(2) for Zi , very close to the expected 
value, and a value x = 0.274(9) for the exponent. Using 
now equation (J6j> and the estimate for the A t from fig- 
ure H3] we obtain an estimate for T t = 1.025(10). Thus, 
the overall estimate for the tricritical point is (A t ,T t ) — 
(1.025(10), 2.925(8)) and is marked with a black rhom- 
bus in figure HU We should note here that, although with 
the power-law fitting ((SJ) we get a very nice estimate for 
the value Aq and also an excellent concurrence between 
numerically estimated transition points and the applied 
law, our attempt above aims only at a numerical quan- 
titative approximation for the main part of the diagram 
and not at the correct asymptotic behavior at its ends. In 
fact, is well known that the asymptotic approach of phase 
boundaries to the sections with the axis follows specific 
power-law behaviors with exponents related to the criti- 
cal exponents describing the transitions in these part of 
the diagram [3Tj . 



4 Conclusions 

In the present paper we have performed a careful numer- 
ical investigation of the critical properties of the Blume- 
Capel model embedded in the triangular lattice. By ap- 
plying an extensive two-stage Wang-Landau entropic sam- 
pling we have studied the finite-size scaling properties of 
the model for several values of the crystal field in both 
its first- and second-order phase transition regime and 



we have estimated with high accuracy transition temper- 
atures and critical exponents. 

In particular, for the regime of second-order phase tran- 
sitions we have verified the theoretical expectation that 
the Blume-Capel model belongs to the universality class 
of the simple d — 2 Ising model, whereas for the regime of 
first-order phase transitions we have identified the most 
characteristic features signaling a first-order phase tran- 
sition and we have estimated the surface tension and la- 
tent heat of the transition using the method of Lee and 
Kosterlitz. Finally, using our estimates for the transition 
temperatures and a novel scheme that takes advantage of 
the logarithmic scaling behavior of the specific heat, we 
have proposed an approximation of the phase diagram of 
the model in the overall (A, T)-plane and we have refined 
previous mean-field-type estimates for the tricritical point. 

Closing, we would like to leave as an open research 
challenge the investigation of the phase diagram of the tri- 
angular Blume-Capel model via a multi-parametric Wang- 
Landau method that would give direct access to the esti- 
mation of the phase diagram. Such an analysis would be 
complementary to the present contribution and has been 
already successfully performed by Silva et al. [IU] for the 
case of the square lattice Blume-Capel model. 

The author would like to thank Professor A.N. Berker for useful 
discussions and a critical reading of the manuscript. 
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